%This is a function that draws a variance matrix from inverted wishart
%distribution. I follow the method proposed by Rossi et. al.(2005) p.43-44
function y=iwhish(sigma, df)

%sigma is the variance parameter
%df is the degree of freedom

%K is the dimension of the variance matrix, also the number of parameters
K=length(sigma);


%step 1, get the upper cholesky roots of the inverse of sigma

U=chol(inv(sigma));

%step2, simulate from the standard wishart to get T

T=zeros(K,K);

for i=1:K
    for j=1:K
        if i==j
            T(j,i)=(chi2rnd(df+1-j))^(.5);
        end
        if j>i
            T(j,i)=randn;
        end
    end
end

%step 3
C=T'*U;

y=inv(C)*inv(C)';

